
AIAA 2002-3286 

Adjoint-Based, Three-Dimensional 
Error Prediction and Grid Adaptation 

Michael A. Park 

NASA Langley Research Center, Hampton, VA 23681 


32nd Fluid Dynamics Conference 
24-27 June 2002 
St. Louis, MO 


For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics 
1801 Alexander Bell Drive, Suite 500, Reston, VA 20191-4344 



Adjoint— Based, Three-Dimensional Error 
Prediction and Grid Adaptation 

Michael A. Park * 

NASA Langley Research Center, Hampton, VA 23681 


Abstract 

Engineering computational fluid dynamics (CFD) 
analysis and design applications focus on output func- 
tions (e.g., lift, drag). Errors in these output functions 
are generally unknown and conservatively accurate so- 
lutions may be computed. Computable error estimates 
can offer the possibility to minimize computational 
work for a prescribed error tolerance. Such an esti- 
mate can be computed by solving the flow equations 
and the linear adjoint problem for the functional of 
interest. The computational mesh can be modified 
to minimize the uncertainty of a computed error es- 
timate. This robust mesh-adaptation procedure au- 
tomatically terminates when the simulation is within 
a user specified error tolerance. This procedure for 
estimating and adapting to error in a functional is 
demonstrated for three-dimensional Euler problems. 
An adaptive mesh procedure that links to a Computer 
Aided Design (CAD) surface representation is demon- 
strated for wing, wing-body, and extruded high lift 
airfoil configurations. The error estimation and adap- 
tation procedure yielded corrected functions that are 
as accurate as functions calculated on uniformly re- 
fined grids with ten times as many grid points. 

Introduction 

Engineering problems commonly require computa- 
tional fluid dynamics (CFD) solutions with functional 
outputs of specified accuracy. The computational re- 
sources available for these solutions are often limited 
and errors in solutions and outputs are often unknown. 
CFD solutions may be computed with an unnecessarily 
large number of grid points (and associated high cost) 
to ensure that the outputs are computed to within a 
required accuracy. A method to estimate the error 
present in a computed functional offers the possibility 
to avoid the use of overly refined grids to guarantee 
accuracy. 

Unstructured grid technology promises easier initial 
grid generation for novel complex three-dimensional 
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(3D) configurations compared with structured grid 
techniques. The use of unstructured grid technology 
for CFD simulations allows more freedom in adapt- 
ing the discretization of the meshes to improve the 
fidelity of the simulation. Many previous efforts have 
attempted to tailor the discretizations of unstructured 
meshes to increase solution accuracy while reducing 
computational cost. 1-8 

Most of these adaptive methods focus on modifying 
discretizations to reduce local equation errors. These 
local errors are not guaranteed to directly impact er- 
ror in global output functions. These methods, often 
referred to as feature-based adaptation, focus on re- 
solving discontinuities or strong gradients in the flow 
field. Flow features (e.g., shocks) can be in the in- 
correct location due to errors elsewhere in the flow 
field. Also, resolving the flow in a location may have a 
minimal effect upon the output function (e.g., a down- 
stream shock). 

If the flow equations are linearized about the flow so- 
lution, the solution of a linear dual problem can yield a 
direct measure of the impact of local primal (flow equa- 
tion) residual on a selected functional output. The 
combination of the primal and dual problems can also 
be applied to yield a correction to a specified func- 
tional on a given mesh. If a specified error tolerance 
in an output function is required, the cost of comput- 
ing a CFD solution can be minimized by adapting the 
discretization of the problem to directly minimize un- 
certainties in the corrected output function of interest. 
Also, the entire adaptive simulation can be terminated 
when the predicted error is equal to a specified toler- 
ance, preventing the waste of computational work on 
an overly large mesh. 

There are many examples of these techniques in the 
finite element communities. 9,10 Pierce and Giles 11 ap- 
plied these methods to finite- volume discretizations. 
Venditti and Darmofal 12,13 have demonstrated these 
methods for compressible two-dimensional (2D) invis- 
cid and viscous flow solutions. The present study is 
essentially a 3D extension of the methods of Venditti 
and Darmofal 13 that adapts the mesh to reduce uncer- 
tainty in an error correction. Muller and Giles 14 have 
also presented results for a similar approach utilizing 
a different adaptation parameter that directly targets 
the error correction. 

The use of adjoint variables (solution of the dual 
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problem) is an efficient method for computing deriva- 
tives of a functional of interest for gradient-based 
design methods. Some examples of discrete adjoint de- 
sign methods are given in Nielsen and Anderson. 15,16 
The discrete dual problem for adjoint variables can be 
expensive. However, the adjoint solution is already 
available for use during the aerodynamic design pro- 
cess, so it could be employed for simultaneous design, 
error prediction, and grid adaptation. 

The combination of adjoint-based grid adaptation 
and design techniques can yield an attractive tool 
for the aerodynamic design of new configurations. 
Adjoint-based error prediction and adaptation can 
yield smaller meshes than traditional feature-based 
schemes with computable error estimates on out- 
put functions. Design processes require analysis and 
derivative evaluation tools that operate with mini- 
mal human interaction. Robust, automatic adaptation 
techniques enable the increased use of nonlinear flow 
calculations in larger multidisciplinary design frame- 
works. These new techniques will enable efficient 
analysis for existing configurations and expanded ex- 
ploration of design spaces for new configurations. 

Flow Equations 

The FUN3D 17-19 '!' (Fully Unstructured Navier- 
Stokes Three-Dimensional) suite of codes is employed 
in this study. The compressible flow solver employs 
an unstructured finite- volume tetrahedral method for 
conserved variables, Q, i.e., 

Q = [p pu pv pw E] T (1) 

where p is density, u, v, and w are velocity, and E is 
total energy per unit volume. The incompressible flow 
solver employs the following state variables: 

Q = \p u v iv] T (2) 

where p is pressure. The node-based variables Q are 
computed by driving the flow equation residual R to 
steady-state with an implicit point-iterative method. 
The code is able to solve incompressible, Euler, and 
Reynolds-Averaged Navier-Stokes (RANS) flow equa- 
tions loosely coupled to the Spalart-Allmaras 20 one- 
equation turbulence model. The present study em- 
ploys only the Euler and incompressible equations. 
The solution of Q allows the calculation of integral 
quantities / (e.g., lift, drag). To speed execution, the 
problem domain is decomposed and the flow and the 
adjoint problems are solved with a. parallel execution 
scheme utilizing the Message Passing Interface (MPI) 
standard. The FUN3D suite of codes is being extended 
to the HEFSS (High Energy Flow Solver Synthesis) 21 
modular framework of FORTRAN 90 shared libraries. 

thttp: //fun3d. larc .nasa.gov 


Adjoint Equations 


After the flow solution is known, the discrete ad- 
joint equations 15, 18,19 are solved to complete the dual 
problem. The first step is to linearize the flow equation 
residual R and output function / with respect to the 
flow solution Q. After this linearization, an adjoint 
variable A is solved for each of the flow equations. 

An abbreviated derivation, adapted from Taylor et 
al., 22 is below. The chain rule for the linearized output 
function is 


df _ ( df \ T dR 
3Q V<9rJ 9Q 


(3) 


The adjoint variable A is defined as the effect of the 
flow residual on the output function: 


2i 

dR, 


= A 


A set of linear equations is solved to find A: 


0R\ 

dQJ 


T 


A = 



(4) 

(5) 


After the flow solution is known, this set of linear equa- 
tions is solved with GMRES. 23 See Refs. 15,16, and 19 
for details. A implicit point-iterative time-marching 
method is employed to compute the adjoint solution 
for the high lift, configuration. 24-26 


Error Correction 

The error prediction and correction scheme is taken 
from Ref. 13. With a solution on a mesh of reasonable 
size Q°. it is desirable to predict the value of an out- 
put. function evaluated with a solution on a much finer 
mesh /( Q*). This prediction can be computed with- 
out the solution on this finer mesh when the adjoint 
solution on this reasonable mesh A° is utilized. The 
full derivation of the error correction term is available 
in Refs. 11, 13. 14. An abbreviated derivation is pre- 
sented by expanding the Taylor series about /( Q°), 
i.e., 


/( Q*) = /( Q°) + 



R(Q°)) + ... 


( 6 ) 

Employing eq. (4) and assuming that the residual 
on the much finer mesh is zero yields an improved es- 
timate for the functional of interest f est : 


2i 

OR 


o 


A° 


(7) 


R(Q*) = 0 (8) 

/( Q*) « fest = /( Q°) - (A°) t R(Q°) (9) 

To improve the prediction of the functional f est , Q" 
and A° can be interpolated to an embedded mesh. In- 
terpolation is performed in two ways for this study: 
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a linear interpolation (Q L ,A L ) and a higher order in- 
terpolation (Q h ,X h ). Details of this interpolation and 
the construction of this embedded mesh are in the “In- 
terpolation Techniques” and “Embedded Mesh” sec- 
tions. Substituting these interpolated quantities into 
eq. (9) yields the linear and higher order functional 
estimates / e L st and f^ gt : 

fest '**fi Q L ) - (A l ) t R(Q l ) (10) 

feat = /( Q H ) ~ (X H ) T R(Q H ) (11) 

Adaptation Parameter 

The adaptation parameter, also from Ref. 13, is in- 
tended to specify a grid spacing modification to reduce 
the uncertainty in the computed error prediction. The 
grid is not modified to directly reduce the computed 
error prediction (as in Ref. 14) because an estimate for 
the functional including this error term can be com- 
puted with eq. (9). Instead, targeting the uncertainty 
in this computed quantity is more effective and im- 
proves the robustness of the adaptive process. The 
error correction (eq. (9)) including the uncertainty in 
the dual solution is 


/( Q°) - /( Q*) « (A°) t R(Q°) + (A* - A°) r R(Q°) 

( 12 ) 

The uncertainty in the computed error correction is 

feat ~ /(Q*) « (A* - A°) T R(Q°) (13) 

The relation of the primal and dual problems 11,13 
yields another expression for the error correction un- 
certainty 

(A* - A°) T R(Q°) = R a (A°)(Q* - Q°) t (14) 


Where Ra(A) is the residual of the dual problem: 


Ra(A) 



(15) 


A computable term is found by using the interpolation 
error of A to replace (A* — A 0 ) and the interpolation 
error of Q to replace (Q* — Q°). The higher order 
interpolate for Q 11 and A° is employed to improve pre- 
diction in place of the linear interpolate in Ref. 13. 
The interpolation error is expressed as the difference 
in the high-order and linear interpolated values: 


(A* 

(Q* 


A 0 ) « ( X H - X L ) 

(16) 

Q°)«(Q ff -Q i ) 

(17) 


The average of the absolute values of the two uncer- 
tainty terms in eq. (14) yields the adaptation intensity 
I, which is computed for each equation on each em- 
bedded node: 


(A ff - A L ) T R(Q H )| + |R A (A ff )(Q H - Q L ) T | 

2 

(18) 


The intensity I is therefore the average of a primal 
residual weighted with a dual solution interpolation 
error and a dual problem residual weighted with a 
primal solution interpolation error. This form of the 
adaptation intensity tends to focus on the nonlinear 
contributions to the function error, which increases ro- 
bustness of the adaptation method. 

Error Correction and Adaptation 
Process 

The error correction and adaptation process begins 
with an initial tetrahedral mesh, which can come from 
any mesh generation system. The state variables are 
computed as the nonlinear solution to the flow equa- 
tions on the initial mesh. The adjoint variables are 
then computed with the linearized flow equations at 
the flow solution. These flow and adjoint solution pro- 
cedures employ a parallel execution scheme. Then the 
global problem domain is reconstructed to facilitate 
the creation of a finer, embedded grid with interpo- 
lated primal and dual solutions. 

Embedded Mesh 

To compute the error prediction and the adapta- 
tion parameter a globally embedded or h-refined mesh 
is created. To construct the embedded mesh, a new 
node (open circle) is inserted at the midpoint of each 
existing edge that connects existing nodes (closed cir- 
cles); see Fig. 1(a) and 1(b). Each existing tetrahedron 
is subdivided to reconnect these new nodes with eight 
interior tetrahedra. (Each of the existing boundary 
faces is also divided into four triangles.) 



a) Original tetrahedron. 



b) Embedded tetrahe- 
dron. 


Fig. 1 Tetrahedron embedding process. 


The four new tetrahedra constructed at the corners 
of the existing tetrahedron have the same shape as 
the original tetrahedron but are smaller in size. The 
construction of the four corner tetrahedra leaves an 
interior volume with eight faces, which is subdivided 
into four tetrahedra. The four interior tetrahedra have 
three unique configurations. The configuration with 
the lowest maximum dihedral angle is selected. 

The new nodes are placed at the midpoints of edges 
during the mesh embedding process. The embedded 
nodes on the boundaries of the mesh may no longer 
remain on the original surface definition of the model. 
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When the grid is adapted to improve the discretiza- 
tion, the surface fidelity of the mesh is maintained with 
boundary node projection. 

Interpolation Techniques 

The primal and dual solution variables y, which are 
all elements of Q and A, are now interpolated to this 
embedded mesh. The value of the solution at each 
of the existing nodes is directly copied into the corre- 
sponding nodes of the embedded mesh. Each of the 
solution variables is interpolated in two ways to form 
the linear and the higher order reconstruction for the 
new nodes. The higher order reconstruction of the so- 
lution for the new nodes requires the computation of 
least-squares gradients 17 at the existing nodes using 
the existing mesh and solution. To simplify the 3D 
interpolation, the interpolation is performed by inde- 
pendently examining each existing edge in the original 
mesh; thus the interpolation problem becomes one- 
dimensional along each existing edge. Each existing 
edge has an existing node at each edge endpoint and 
a new node at the edge midpoint (see Fig. 2). 



Fig. 2 New node and existing edge. 

An edge has two 3D endpoints: Xi and X 2 . The 
vector that represents the length and direction of the 
edge is Ax = X 2 — xi. The 3D least-squares gradient of 
the solution Ay can be projected to a total derivative 
along the edge to facilitate interpolation by 

y = S7y T Ax (19) 

The new node interpolation y 3 / 2 can be expressed 
as a combination of the solution values and the deriva- 
tives at each endpoint (y 3 , yi, y 2 , and y 2 ). The linear 
interpolation y 3 / 2 is the average of the two end nodes 


Error Correction and Adaptation Parameter 

Once the new embedded grid is constructed with 
Q L , A L , Q h , and X H , it is partitioned to allow paral- 
lel calculation of the functional and flow and adjoint 
equation residuals. The flow and the adjoint equations 
are not iterated or solved on this embedded grid; the 
flow state and adjoint variables are interpolated from 
the original mesh. Therefore, the only computational 
costs on this larger embedded grid are function evalu- 
ations, flow and adjoint residual evaluations, and dot 
products of vectors. The linear and higher order er- 
ror correction term, eq. (10) and (11), is computed at 
each node on the embedded mesh and summed over 
the entire mesh for all flow equations. 

The adaptation intensity, eq. (18), is also computed 
at each embedded mesh node. At each node that is 
also present in the original mesh, the computed inten- 
sity is zero due to the chosen interpolation technique. 
The values of the lower and higher order interpolation 
schemes are the same at these existing nodes; thus 
( Q h — Q l ) and ( X H — X L ) are exactly zero. 

To specify the grid adaptation on the original mesh, 
the adaptation intensities must be reduced from the 
embedded mesh to the original mesh (Io). The new 
nodes on the embedded mesh all he on existing edges 
of the original mesh (see Fig. 1(b)). Therefore, to con- 
struct Io the original mesh is examined one edge at 
a time (see Fig. 2). One half of the intensity com- 
puted at each new node (which is at the midpoint of 
these original edges) is added to each existing node at 
the endpoints of these edges. The intensities are also 
summed over the equations at this point, resulting in 
one intensity value for each original node. 

The adaptation parameter, which has been reduced 
to the original mesh, is summed to find the global 
intensity I g = YI Io • The number of nodes in the orig- 
inal mesh n and the user-specified error tolerance t are 
combined to scale the adaptation intensity; that is 

h = ( 22 ) 


2/3/2 


2/1 + 2/2 
2 


( 20 ) 


The higher order interpolation y 3 / 2 is found with a 
cubic fit of the endpoint data that is evaluated at the 
midpoint 


H 2/i + 2/2 , 2/i - 2/2 

2/3/2 ~ ^ 


( 21 ) 


Equation (21) is equivalent to a least-squares quadratic 
fit of the endpoint data, that is evaluated at the edge 
midpoint. 

At the completion of the grid embedding and inter- 
polation step, the linear and higher order interpolated 
solutions to the primal and dual problems Q L , A L , 
Q ff , and X H are available to compute the error cor- 
rection and adaptation parameter. 


To perform the grid adaptation, the mesh is locally 
enriched in the location of nodes where the scaled 
intensity I s is greater than a value, e.i., unity. Ref- 
erence 13 demonstrates a way to specify a new el- 
ement spacing function. The adaptation procedure 
self-terminates as all elements of I s become less than 
unity (i.e. no nodes are flagged for adaptation). 

Adaptation Mechanics 

The adaptation mechanics utilize three independent 
modules. The first module inserts new nodes into the 
existing mesh and locally reconnects tetrahedra and 
boundary faces to maintain a valid tessellation. The 
second module employs face and edge swapping to im- 
prove the mesh quality. The final module performs 
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grid smoothing and boundary node projection opera- 
tions. 

Node Insertion 

The node insertion method is currently one level 
of selective h-refinement. To start the refinement, all 
the edges surrounding nodes on the original mesh that 
have a scaled intensity I s greater than unity are flagged 
for h-refinement. The set of flagged edges is exam- 
ined tetrahedron by tetrahedron and additional edges 
are flagged in an attempt to maintain grid quality 
(i.e., low maximum dihedral angles, few high-degree 
nodes) . The final set of flagged edges results in tetra- 
hedra with one edge, three edges on one face, or all six 
edges flagged. A tetrahedron with all six edges flagged 
is illustrated in Fig. 1(b). The mesh is then h- refined 
by inserting new nodes on the midpoint of the flagged 
edges and reconnecting these nodes into new tetrahe- 
dra and boundary faces. 

Face and Edge Swapping 

The current post-adaptation grid-improvement 
scheme employs face and edge swapping. 27 The swap- 
ping algorithm minimizes a shape (cost) function (e.g., 
aspect ratio or dihedral angle). This study sought to 
reduce only the cell aspect ratio AR 

1 tetrahedral circumsphere radius 

i4.li — ~ i i i ■ i ~ ( 23 ) 

3 tetrahedral in-sphere radius 

Reconnections of tetrahedra with undesirable shape 
measures are investigated and new local tetrahedra 
configurations with more desirable shape measures 
are selected. Edges on boundary faces can also be 
swapped. To simplify and speed up the edge swapping 
routine, the boundary face information is discarded 
and reconstructed at the end of the swapping process. 
Smart-Laplacian smoothing 28 is used on the interior 
nodes. The actual locations of the boundary nodes is 
not modified in this module; that modification is per- 
formed by the grid smoothing and projection module. 

Grid Smoothing and Projection 

The inserted boundary nodes may not be located 
on the surface geometry of the model to be simulated 
since they were inserted at the midpoints of existing 
edges. A CAD model is employed to describe the 
actual model surface. To regain the surface fidelity 
of the mesh, the newly inserted boundary nodes are 
projected to the model surface with a CAD interface 
package CAPRI (Computational Analysis PRogram- 
ing Interface). 29,30 The projection of these new nodes 
to their location on the CAD surface can result in in- 
verted, invalid tetrahedral elements. 

A grid node-smoothing algorithm is employed to 
facilitate boundary projection without generating in- 
valid elements and to improve the overall quality 
(shape measure) of the mesh. This post-adaptation 
smoothing is similar to Freitag 28 and was initially im- 


plemented by Brasher and Park. 31 This initial imple- 
mentation has been extended and coupled to CAPRI, 
which utilizes native CAD point projection routines. 
Nodes on the boundary are smoothed by moving the 
nodes in CAD ( u,v ) parametric space to improve the 
shape measure of adjacent tetrahedra. 

As the nodes are projected, the neighboring tetra- 
hedra are tested for validity. If invalid tetrahedra 
resulted from the projection, the projection distance 
of the boundary nodes is reduced until the neighboring 
tetrahedra are valid. Then the nodes in the neighbor- 
hood of the projected node are smoothed to improve 
a quality measure of the adjacent tetrahedra. The 
boundary points are then moved into the fully pro- 
jected position in a number of iterative cycles. 

It is anticipated that grid smoothing in the neigh- 
borhood of projected nodes may not adequately re- 
gain surface fidelity of anisotropic meshes. A grid- 
movement scheme may be required as in Ref. 16. An- 
other possibility is a 3D version of mesh restructuring 
as in Ref. 32. 

Adaptation Module Interaction 

The current selective h-refinement technique often 
creates high-degree nodes on the border of the adapted 
regions. The smoothing algorithm is currently unable 
to improve elements that are adjacent to high-degree 
nodes. The edge and face swapping techniques effec- 
tively improve shape measures and reduce the number 
of high-degree nodes, facilitating projection and node 
smoothing. 

These three adaptation modules where developed 
independently to facilitate a quick initial implementa- 
tion and to investigate the strengths and weaknesses 
of each technique. Merging the abilities of these three 
separate modules will allow for more flexible modifica- 
tions of grids (e.g., point insertion, point removal, and 
anisotropic elements). 4,6 

Results 

Adaptation results are shown for a wing, wing-body, 
and high lift configurations. The wing is simulated 
with incompressible and transonic flow conditions. 
The wing-body and high lift configurations are sim- 
ulated with subsonic flow. 

Initial Mesh Generator 

The initial meshes for these error prediction and 
adaptation studies are generated with the FELISA 
mesher 4 connected to CAD geometry by CAPRI 
through the GridEx 33 framework. FELISA is a Delau- 
nay mesh generator with an advancing-front method 
for inserting nodes. The GridEx framework is cur- 
rently being developed at NASA Langley Research 
Center to link various grid generation and adaptation 
strategies to geometry through CAPRI. This frame- 
work is also utilized in a batch mode to perform uni- 
form grid refinement studies. 
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a) Initial ONERA M6 mesh. 



b) ONERA M6 mesh adapted to incompressible drag. 


Fig. 3 Initial and adapted ONERA M6 meshes. 

Drag Adaptation — Incompressible ONERA M6 

The initial mesh for an ONERA (Office National 
d’Etudes et de Recherches Aerospatiales) M6 wing 
with 5227 nodes is shown in Fig. 3(a). The mesh 
has extremely coarse spacing, especially at the trailing 
edge and is intended to resolve the surface curvature of 
the leading edge and the wingt.ip and not any particu- 
lar flow features. The spacing function for this mesh is 
specified manually and is intended to be representative 
of an automated curvature or maximum chord-height 
specification. The CAD geometry is represented with 
the CAPRI native kernel FELISA with a part gener- 
ated from a surface ICES definition. 

The initial ONERA M6 mesh was used in the grid 
adaptation process with incompressible flow at an an- 



1000 10000 100000 1e+06 

Number of Nodes 

Fig. 4 Coefficient of drag for the ONERA M6 
adapted to incompressible drag. 

gle of attack of 0 deg. Directly computed drag and 
estimates of drag for an ONERA M6 wing as a func- 
tion of number of nodes is shown in the Fig. 4 log-log 
plot. The adaptation and error correction results are 
shown for a drag error tolerance of 0.001. The directly 
computed drag on the adapted meshes is represented 
by the solid lines with circular symbols. The error- 
corrected drag calculated with the linear interpolated 
solution f 3St is represent ed by a dashed line and square 
symbols. The estimated functional calculated with the 
higher order interpolated solution f^ st is represented 
by a dashed line and diamond symbols. The correct 
drag is zero because of non-lifting, subcritical, invis- 
cid flow. Therefore, the y-axis denoted “Coefficient of 
Drag” is also the error in drag. 

The triangle labeled “Second-Order Slope” in Fig. 4 
illustrates second-order spatial convergence. The 
adaptive-grid method results in drag calculations 
that converge at a much higher rate than the 
asymptotic convergence rate of uniform refinement 
(second-order). The adaptive procedure correctly self- 
terminated when the drag error of the adapted flow 
solution reached the user specified error tolerance (dot- 
dash line). 

The final grid (454 thousand nodes) after five cycles 
of grid adaptation to incompressible drag is shown in 
Fig. 3(b). The adaptation process clustered grid points 
at the leading and trailing edges of the wing. Points 
are also clustered in the neighborhood of the stagna- 
tion stream line. 

The tetrahedra shape measure AR (eq. 23) is min- 
imized by the mesh improvement algorithm. The 
boundary node smoothing algorithm is intended to op- 
timize the shape measures of the tetrahedral elements. 
Therefore, the shape measures of the boundary faces 
depicted in this surface plot may not be optimal. 
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b) Richardson extrapolation. 

Fig. 5 Coefficient of drag for the uniformly refined 
ONERA M6 at 0.84 Mach. 


Drag Adaptation — Transonic ONERA M6 

A uniform refinement of the ONER A M6 wing mesh 
is computed at 0.84 Mach and an angle of attack of 3 
deg. The drag directly computed by the flow solver 
is shown with linear and higher order interpolated er- 
ror corrections as a function of the number of nodes 
in Fig. 5(a). The extrapolated (grid-converged) drag 
value for Fig. 5(a) was estimated with Richardson ex- 
trapolation from Fig. 5(b). These meshes have the 
same spacing function as Fig. 3(a) globally modified 
with a scalar to uniformly reduce the element spac- 
ing. These grids were generated with the batch version 
of GridEx using the FELISA mesher and CAPRI for 
CAD geometry access. 

Figure 5(b) shows drag and estimates of drag as 
a function of element size for the uniform grid re- 

7 




Fig. 6 Coefficient of drag for the drag adapted 
ONERA M6 at 0.84 Mach. 


finement. of the ONERA M6 wing. A representative 
element length h was estimated by computing the cube 
root of the number of nodes. This length was nor- 
malized by the estimated length of the 624 thousand 
node mesh ho- The symbols are drag computed by 
the flow solver and error corrected values. A linear 
fit of the data at (ho/h ) 2 = 1.0 and (ho/h ) 2 = 1.7 
is used to estimate the grid-converged answer for all 
three schemes. All three schemes indicate a similar 
grid-resolved value. An additional flow solution (1.2 
million nodes, (ho/h ) 2 = 0.6) is shown to verify the 
accuracy of the linear fit; it was not used to construct 
the linear fit of the computed drag. An improved in- 
terpolation scheme for the error correction would yield 
a superconvergent functional estimate as in Ref. 11. 

The initial coarse mesh shown in Fig. 3(a) is em- 
it 
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Fig. 7 Initial ONERA M6 upper wing surface. 


ployed in the error prediction and grid adaptation 
procedure with two different adaptation methods. The 
coefficient of drag is plotted as a function of mesh size 
in Fig. 6. The user specified error tolerance in drag 
is 0.0019. The uniformly refined flow solution from 
Fig. 5(a) is shown with the adapted grid flow solu- 
tion and higher order error prediction of the adapted 
grid in Fig. 6. Figure 6(a) demonstrates a single level 
of h-refinement, for all nodes with a scaled adaptation 
intensity I s greater than one. Figure 6(b) shows h- 
refinement for I s greater than one and a recursive call 
to the adaptation mechanics for I s greater than 75, 
yielding two levels of h-refinement at each adaptation 
cycle. The initial convergence of the function is better 
in Fig. 6(b) than Fig. 6(a). This improvement in func- 



a) Upper wing surface mesh. 



Fig. 8 Final ONERA M6 upper wing surface, 
adapted with one level of h-refinement. 


tion convergence is illustrating the limitations of using 
a single level of selective h-refinement as the adap- 
tive node-insertion procedure. The use of two levels of 
h-refinement better approximates a continuous vari- 
ation in element size. The two adaptation methods 
converged to similar meshes and drag values. 

The upper wing surface grid and Mach contours 
of the initial flow field computed on the mesh from 
Fig. 3(a) is shown in Fig. 7. The shocks in this 
initial grid are poorly resolved. The mesh (357 thou- 
sand nodes) and Mach contours of the ONERA M6 
adapted to drag with one level of selective h-refinement 
is shown in Fig. 8. The mesh (374 thousand nodes) 
and Mach contours of the ONERA M6 adapted to 
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a) Upper wing surface mesh. 



Fig. 9 Final ONERA M6 upper wing surface, 
adapted with two levels of h-refinement. 


drag with two levels of selective h-refinement is shown 
in Fig. 8. The final meshes and solutions are simi- 
lar for both of the adaptation methods. The adaptive 
procedure strongly clustered nodes at the leading and 
trailing edges of the wing and lightly clustered nodes 
at the shock location on the upper surface. Feature- 
based adaptations of this configuration in Ref. 7 and 8 
focused on the leading edge and shock locations, but 
not the trailing edge. 

Lift Adaptation — EET 

The EET (Energy Efficient Transport.) 34 initial 
coarse mesh is shown in Fig. 10. The initial grid spac- 
ing distribution is specified manually to resolve the 




Fig. 11 Coefficient of lift for the lift-adapted EET 
at 0.40 Mach. 

surface details of the fuselage, wing leading edge, blunt 
trailing edge, and the wingtip. The geometry is repre- 
sented with a Parasolid CAD kernel accessed though 
the CAPRI application program interface (API). 

The initial coarse mesh shown in Fig. 10 is employed 
in the lift error prediction and grid adaptation proce- 
dure at 0.40 Mach an angle of attack of 2 deg. The 
lift coefficient is plotted as a function of mesh size 
in Fig. 11. The adaptive procedure has a user speci- 
fied error tolerance of 0.1 for the lift coefficient error. 
The uniformly refined lift calculation is shown with the 
adapted grid lift calculation and error predictions on 
the adapted grid in Fig. 11. The Richardson extrapo- 
lation value is not shown because a reasonable linear 
fit of the last three points was not possible. The lift 
coefficient is calculated on an adapted mesh one-tenth 
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5 



Fig. 12 Initial 30P-30N grid. 


a) Initial 30P-30N grid. b) Lift-adapted 30P-30N 

grid. 

Fig. 13 Original and adapted 30P-30N symmetry 
plane grids. 

the size of the uniformly refined grid. The adapted lift 
error is well within the user specified tolerance. 

Lift Adaptation — 30P-30N Airfoil 

The McDonnell Douglas 30P-30N airfoil initial 
coarse mesh is shown in Fig. 12. The 30P-30N airfoil 
is extruded between two symmetry planes. The near- 
plane has been removed to improve visualization. The 
geometry is represented with a Parasolid CAD kernel 
accessed though the CAPRI API. This configuration 
is the subject of a recent 3D CFD study. 35 

The geometry and initial coarse mesh (113 thousand 
nodes) shown in Fig. 12 is employed in the error pre- 
diction and grid adaptation procedure at 0.20 Mach 
and an angle of attack of 16.3 deg. The lift adaptive 
procedure has a user specified error in lift, of 0.25. The 
uniformly refined grid flow solution, adapted grid flow 
solution, linear error prediction, and higher order er- 
ror prediction are shown in Fig. 14. The extrapolated 
coefficient of lift value was computed with a Richard- 
son extrapolation of the finest two uniformly refined 
solutions. The original symmetry plane grid and the 
symmetry plane grid (832 thousand nodes) after two 



100000 1e+06 

Number of Nodes 

Fig. 14 Coefficient of lift for the lift-adapted SOP- 
SON at 0.20 Mach. 

cycles of adaptation is shown in Fig. 13. 

Table 1 shows the aspect ratio AR (eq. 23) for the 
initial and the adapted grids. The AR is the cost 

Table 1 Shape measure for the 30P-30N adapta- 
tion 


Cycle 

Aspect ratio AR 

Face angle 


max 

max 

1 

8.1 

154.3 

2 

5.0 

158.1 

3 

7.5 

167.8 


function for the grid-improvement optimizer. The face 
(dihedral) angle is not directly controlled but could be 
added as a constraint. 

Conclusion 

The initial implementation of an adjoint-based error 
correction and adaptation method has been demon- 
strated in three dimensions. With a given flow and 
adjoint solution, the error correction for a functional 
and adaptation intensity term have been described. 
The adaptation intensity was formulated to reduce the 
uncertainty in the error correction of a global func- 
tional and not a local error estimate as in a feature- 
based adaptation scheme. The adaptive procedure 
automatically terminates the simulation when an user 
specified tolerance is satisfied. The error remaining in 
the simulation at termination was always within the 
user specified tolerance, although sometimes the sim- 
ulation was overly accurate. 

A wing configuration was adapted to reduce drag 
error in incompressible and transonic flow. The drag 
computed by this adaptation and error correction 
method was shown to be as accurate as direct flow 
calculations using larger uniformly refined grids. The 
initial convergence of adaptation procedure improved 
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with two levels of h-refinement at each adaptation 
cycle. Lift adaptations of Parasolid CAD models of 
wing-body and high lift configurations demonstrate 
the utility of this adaptive methodology on complex 
geometries. The lift of the wing-body configuration 
was computed on an adapted grid that is one-tenth 
the size of an uniformly refined grid. 
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